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Abstract 

In this paper we provide the derivation of a super compact pairwise model with 
only 4 equations in the context of describing susceptible-infected-susceptible (SIS) 
epidemic dynamics on heterogenous networks. The super compact model is based on 
a new closure relation that involves not only the average degree but also the second and 
third moments of the degree distribution. Its derivation uses an a priori approximation 
of the degree distribution of susceptible nodes in terms of the degree distribution of 
the network. The new closure gives excellent agreement with heterogeneous pairwise 
models that contain significantly more differential equations. 
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1 Introduction 


While networks have provided a new modelling paradigm for population dynamics mmm, 
these are still used in conjunction with mean-field models of various types. The most 
frequently used and well-known mean-field models for network epidemics are the degree- 
based mean-field (DBMF) model, also known as heterogeneous mean-field [3] [3] and pair¬ 
wise model Eiimtzi- Both continue to provide a productive framework for approximating 
expected values of random variables emerging from explicit network-based stochastic sim¬ 
ulations in different contexts and networks with different properites. The major advantage 
of such mean-field models stems from the fact that often these allow us to analytically 
determine quantities such as the basic reproduction number, final epidemic size or en¬ 
demic equilibrium H El- Such analytic expressions then lead to a significantly better 
understanding of the interplay between network and disease characteristics. 

Pairwise models have originally been introduced in the context of mathematical ecol¬ 
ogy m followed by natural extensions to epidemiology [B]. The original simple model for 
undirected and unweighted networks has been subsequently extended to networks with het¬ 
erogenous degree [3j, directed networks m, weighted networks m, networks displaying 
motifs [12j and even combined with the edge-based compartmental modelling framework 
for an even more compact treatment [7]. 

The closure in the most basic or fundamental pairwise model is based on the assump¬ 
tion on homogeneity of the degree distribution, i.e., all nodes have approximately the same 
degree n. Hence the traditional pairwise model cannot be applied for graphs with hetero¬ 
geneous degree distribution, such as bimodal graphs or networks with power law degree 
distribution. This is shown in Figured] For heterogeneous networks, a corresponding pair¬ 
wise model was introduced in |9j. This gives excellent agreement with simulations for all 
configuration-like random networks m , see Figure [T] The heterogeneous pairwise model 
consist of order N 2 differential equations, where N denotes the number of nodes in the 
network. An approximation of pairs leads to a simpler system, called compact pair-wise 
model that consist of only order N equations [7] and still gives very good agreement with 
simulations, see Figured] 

The aim of this paper is to introduce an even simpler model with only four equations 
that performs well for large heterogeneous networks. The system is derived from the 
compact pairwise model by introducing a further approximation, and using a closure 
relation that contains not only the average of the network’s degree distribution but also 
its second and third moments. 

2 Derivation of the super compact paiwise (PW) model 

2.1 Pairwise model for homogenous networks 

We start from the exact PW model. For the SIS epidemic on an arbitrary undirected 
network the expected values of [5], [/], [SI], [II] and [SS] satisfy the following system of 
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differential equations 


[S\ = 7 [I\-r[SI\, (1) 

[I] = t[SI] — 7[i], (2) 

[SI] = 7 ({II\-[SI})+t([SSI]-[ISI}-[SI}), (3) 

[SS\ = 2 7 [SI]-2t[SSI], (4) 

[//] = -2 1 [II} + 2t([ISI] + [SI}), (5) 


where [X], [XY] and [XYZ] denote the expected number of nodes in state X, edges in state 
X — Y and triples in state X — Y — Z , respectively, with counting according to all possible 
edge directions. This effectively means that for undirected networks [XY] = [YX], [XX] 
is double the number of unique edges in state X — X, and similarly X — Y — X accounts 
twice for a unique X — Y — X triple, where X, Y E {S, /}. This system is derived directly 
from master equations in jTl] and hence exact. We note that some of the equations can 
be omitted by exploiting conservation identities, such as [S'] + [I] = N. 

It is well known that in order to transform Eqs. (HI) - (EJ) into a self consistent solvable 
system closures need to be applied in order to break dependency on higher order moments. 
Particularly useful are closures at the level of triples. As it is well-known, the simplest 
closure is 

n — 1 [AS] [SI] 

n [S] ’ W 

where n = ( k) is the average degree of the network, and A stands for S or I. This closure 
leads to the traditional pairwise system 



[S\p 

= 7 [I\p~r[SI] p , 

(7) 

\I\v 

= t[SI] p -'y[I] p , 

(8) 

[S/] P 

n |ojp 

(9) 

[SS] P 

_ ^n-l[SI\ p [SS\ p 

~ 2 7[MJ p ^ n [5]p , 

(10) 

\II\v 

n - 1 [SI]l 

= -2 7 [//] p + 2r-^ + 2r[% 

(11) 


Here the subscript p is used to emphasize that the solution of this system is different from 
the exact values of the expected variables. As Figure 1 shows, this system cannot capture 
network heterogeneities, hence closure ([6]) needs improvement. 


2.2 Pairwise models for heterogenous networks: the heterogeneous, pre¬ 
compact and compact pairwise models 

The problem with closure ([6]) is that it assumes that each node has degree n, which is 
obviously a crude approximation for heterogeneous networks. This has led to several 
heterogeneous mean-held models, where the state space is much extended to account for 
the expected number of nodes in different states and with a given degree, e.g., [S'/ C ](t) and 
[Ik] ( t ) for the expected number of susceptible and infected nodes of degree k , respectively. 
These new variables will induce or require further variables at pair level, such as [Skli](t) 
which denotes the expected value of the number of edges connecting susceptible nodes of 
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degree k to infected nodes of degree l. In this spirit, the following heterogeneous models 
were developed in historical order: 

• heterogeneous pairwise model [9], 

• pre-compact pairwise model [9] and 

• compact pairwise model [7]. 

Instead of presenting the systems of differential equations of these models and working 
from the most explicit or complex to the more compact one, we start from the simplest 
model and show in an intuitive way how the more sophisticated models arise. Since closure 
© uses the degree of the middle node, it is useful to express the triple as 

K 

[AS I] = ^[AS k I], 

k =1 


where the different degrees occurring in the graph are k = 1, 2,... , 77. The closure for the 
the triples in the right hand side can be written as 

k-l[AS k ][S k I] 

k [5 fc ] ’ j 

In order to use this closure in the exact system ([I])-© one needs differential equations for 
[5^], for [Ski] and for [5^5]. The exact differential equations for [S k ] are 



[S k ] = l[h ]~r[S k I], A; = 1,2,..., 77, (13) 

where the substitution [/*.] = Nk — [5&] can be used. The simplest heterogeneous model [7] 
uses only [Sk] as new variables and introduces an algebraic expression that approximates 
[Ski] and [SkS] in terms of [Sk], [5/] and [55] as follows: 


[S k I] 


[SI\ 


k[Sk] 

Eh i[Si]’ 


(14) 


which can be interpreted as showing that the ratio of the number of edges connecting 
degree k susceptible nodes to infected nodes and the number of SI edges is almost the 
same as the ratio of the number of stubs starting from degree k susceptible nodes and the 
total number of stubs starting from susceptible nodes. Using this approximation, closure 
m can be simplified as given below 


[AS k I] 


k - 1 [A5fc][5fci] k - 1 [AS][SI]k 2 [S k ] _ [AS][SI]k(k - l)[5 fc ] 


k 


[ 5 fc ] 


k 


SI 


SI 


(15) 


N 

where 5 i = E k[S k ] is the first moment of the distribution of susceptible nodes. This 
k =l 

leads to the so-called compact pairwise model, in which the variables are: [57], [55], 
[II] and [Sk] for k = 1,2,... ,77, i.e., it contains 77 + 3 differential equations. In fact, 
the system consists of equations (Tdl) . and ©-© with the above mentioned closures and 
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approximations, namely <HH> and (1151) . Thus it takes the form 


[4] C 

= 7 [4]c-rfc[Sd c ^, 

(16) 

[4] c 

= rk[S k ] c [ -^-'y[Ik]c, 

dS 

(17) 

[SI]c 

= 7 ([II] C - [SI\c) + r([4S] c - [ST\ c )[Sr\ c P - t[SI] c , 

(18) 

[SS] C 

= 2 7 [5J] c - 2 t[SS] c [SI] c P, 

(19) 

[II\c 

= 2 t[SI] c - 2 7 [//] c + 2 r[SI] 2 c P, 

K K 

= j2 k [Sk] c , p=-^^( k -m s k]c 

(20) 

S s 

(21) 


k =1 s k =1 


Here the subscript c, referring to the word ‘compact’, is used to emphasize that the solution 
of this system is different from the exact expected values. 

The next level of complexity is represented by the pre-compact pairwise model, in 
which the variables [S k I] and [S k S] are kept as independent variables and differential 
equations for these are written down. Thus the systems can be formulated in terms of 
variables such as, [4c]) [4fc<S]> [Ski], [44 and [4-1] > i.e., resulting in a total of 5 K variables. 
This can be done by considering the closure introduced in [9] which is 


[A n Bm] 


[A n B][A n B] [N n N m ] E q qWg] 
[AB] n[N n ]m[N m ] 


( 22 ) 


where N k denotes the number of nodes of degree k. It is wroth noting that this system is 
not able to account for preferential mixing. 

The most complex system, which we call heterogeneous pairwise model, uses all com¬ 
binations of pairs as variables, namely [4t4]> [44] and [4-4 Hence, it consists of 2 K 2 
differential equations. At the price of having a system with the number of equations of 
quadratic order, we do not need any extra approximations (besides the closures), such as 
(1141) in the compact pairwise model, or (1221) for the pre-compact pairwise model. Without 


explicitly including the closures, the most complex system can be written as 

[S k ] = —t ^ l[S k Ii] + 7[4], 

(23) 

[4] 

= + T J2i[s k ii]- 7(4], 

(24) 

[S k Si] 

= -T ^2 ( [ImSkSl] + [SkStlm ]) + 7 ( [Skh] + [44]), 

777 , 

(25) 

[Skh] 

= +t'Y' j ([SkSlIm] - [. ImSkh ])- (t + 'y)[S k Il] +7[44], 

777 . 

(26) 

[44 

= +T Y, ( [ImSkh] + [IkSlIm ]) + T ([S k h] + [44]) ~ 2 7 [4 h], 

m 

(27) 


with all subscripts going from 1, 2,... , K. 

2.3 Super compact pairwise model with heterogeneous triple closure 

We now show that the network heterogeneity can be captured by a small system, containing 
only four differential equations, just as in the simplest pairwise model. Consider a random 
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network with degrees d\, d> 2 , ■ ■ ■, dx and denote the number of nodes of degree dk by 
for k = 1,2,..., K, i.e., iVi + IV 2 + ... + Nx = N. We note that denoting degrees as dk 
instead of k will prove to be advantageous in the derivation below. The degree distribution 
of the graph is then given by pk = ^. The average degree and the second moment of the 
degree distribution are 

K K 

w = d ^ < fe2 > = nT, d l N fc- ( 28 ) 

k =1 k =1 

In order to arrive to our new even more simplified system, the super compact PW 
model, we start from a triple and the closure given in (1151) 

MS/] = 5>SU] « - l)[S t ] = [AS][S/]^h, 

k= 1 k =1 


K 

where we used closures m and m, and where S 2 = ^|[£fc] is the second moment 

k =1 

of the distribution of susceptible nodes. Thus in order to use this closure in the exact 
system (|T |) - (|5|) one needs an algebraic expression of S 2 and S 1 in terms of variables [S'], [/], 
[SI], [II] and [SS] only. Expressing the total number of stubs starting from susceptible 
nodes we get Si = [SI] + [SS] as an exact relation. Thus the problem arises from the 
fact that such an exact relation is not available for the second moment S 2 . Our heuristic 
idea in obtaining a good approximation of is the following. Dividing the equation 

K 

[S] = [<Sfc] by [S] we get that [S/,-]/[S] is a probability distribution. The expected value 
k =1 

of this distribution is known, it is 


K 


12 d k 


[Sk\ 

[S] 


= n s '■= 


[SI] + [SS] 
[S] 


or in other words the average degree of susceptible nodes. Our idea is to use an a priori 
approximating distribution for [S/, : ]/[S] that will be denoted by Sk- This approximating 
distribution satisfies 


si + S 2 + • • • + sk — lj (29) 

disi + d 2 S2 + • ■ ■ + d K SK = ns- (30) 

In order to get an a priori approximating distribution we determined [S^]/[S] numerically 
from the compact pairwise model and compared it to pk = Nk/N, the degree distribution 
of the graph. Numerical results show that these are linearly related, meaning that Sk/pk is 
a linear function of the degree dk- More precisely, Sk/pk can be written as A{t)dk + B(t), 
where A and B are time dependent with this relation assumed to hold for all degrees. This 
allows to deal with the heavily under determined linear system given by Eqs. (EH])-®. 
Introducing the notation qk = Sk/pk the assumption on linearity can be formulated as 

gfc ~ 91 = dK-qi , , 2 K 

dk — d\ d K d\' ’ ' 

This yields an expression for qk in terms of q \, qx and the degrees dk as 

(dx — d\)q k = (dk - d\)qx + (dx ~ dk)qi- 
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Multiplying this equation by p k we get the following relation between s k and p k 

{d K ~ di)s k = p k (d k - di)q K + p k (d K - 4)<?i- (31) 

Observe that q\ and qx can be determined from system (l29lh(l30l) by substituting the 
above expression for s k . Namely, we obtain 

(dx-di) = (ni - di)q K + (d K - n 1 )q 1 , (32) 

(d K -di)n s = (n 2 - nid^qx + (md K - n 2 )qi, (33) 


K 

where n t = ^ d l k p k is the i-tli moment of the degree distribution. (It is more convenient 
k =i 

to use n\ and n 2 instead of (k) and (k 2 ).) Solving the linear system (1321) - (n£Jl) for q\ and 
q K we get 


(n 2 — n\)q\ = n 2 - n\n s + d\{n s - ni), (34) 

(n 2 -n\)q K = n 2 - mns + d K (n s - m). (35) 

Substituting these expressions into (13T1) leads to 

(d K -di)(n 2 -nl)s k = p k {d k -di)(n 2 -nins+dK{ns-ni))+p k {dK-d k ){n 2 -nins+di{ns~ni)). 

Now we are in the position of determining the approximate second moment of the distri¬ 
bution s k . Multiplying the above equation by d 2 and summing from k = 1 to k = K some 
simple algebra yields 


K 

(n 2 - ni)^ d\s k = n 2 (n 2 - n s ni) + n 3 (n s - n). 
k =1 

Note that the third moment n 3 of the degree distribution comes into play. Thus the desired 
quantity S 2 can be approximated as 


K 


S 2 = ^2 d l[Sk] 


k =1 


K 

E 

k= 1 


d 2 [S]s k = [5] 


n 2 (n 2 - ngni) + n 3 (ng - n) 
n 2 - n\ 


Hence using S\ = [SI] + [ SS ] = ng[5] we get 


s 2 -s 1 

Si 


1 _ / n 2 (n 2 - n s ni) + n 3 (ng 

n 2 s [S] V ri 2 — n\ 



Therefore, the new closure relation is 


[AST] 


ns[S] 


f n 2 (n 2 - nsn\) + n 3 (n s - n) 
V n s {n 2 - n\) 



(36) 


We note that in the case of a homogeneous network, where each node has degree n, 
we have ns = n and the average degree is n\ = n. Hence, the expression in the bracket 
simplifies to — — 1. Moreover, the second moment is n 2 = n 2 . Therefore, this term is 

simply (n — 1) and leads to the traditional closure [HST] = • 
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Using the new closure (13611 in the pairwise model (HJ-f© gives the super compact PW 
model in the following form: 


[5], 

= 7[7] s - t[SI] s , 

(37) 

[i] s 

= r[5T] s — 7[/] s , 

(38) 

[s/] s 

= 7 ([//], - [ SI] s ) + r[57] s ([55] fl - [SI] a )Q - r[SI\ a , 

(39) 

[SS] s 

= 2 7 [S/] s - 2 t[S/] s [SS] s Q, 

(40) 

[II]s 

= -2 1 [II] s + 2t[SI] 2 s Q + 2 t[SI] s , 

(41) 


where 

n _ 1 f n 2 (n 2 -ngni) + n 3 (ng-n) \ _ [SI] + [SS] 

nsIS] V n s (n 2 -nf) /’ HS ' [S] 

In the next section we show that this new super compact pairwise model gives an as accu¬ 
rate output as the compact pairwise model, despite of the fact that it contains significantly 
fewer differential equations. 


3 Performance of the new closure for different networks 


As it was shown in the Introduction in Figure (TJ the heterogeneous PW and compact PW 
models give very good agreement with simulations, hence we compare the super compact 
PW model with the new closure to the compact PW model. This comparison will be done 
for different heterogeneous networks. Thus systems (171)- ([TT1) . (fT6l) - (l2Ul) and (1371) - (BT1) will 
be solved numerically and the time dependence of [I] p , [I] c and [I] s are compared, where 
N 

[I] c = [tfc] c is the total number of infected nodes in the compact PW model. The 
k =1 


parameters of the epidemic are fixed at 7 = 1 and r = 3y (k)/(k 2 ). The later is chosen in 
such a way that the ratio of t and its critical value T cr = 7 (k)/(k 2 ) is a given constant. 
Here, this ratio is chosen to be 3, its actual value has only a minor influence on the results, 
generally this need to be greater than 1 to have an epidemic. 

Let us consider first the case of bimodal random graphs, where there are two different 
degrees d\ and d 2 , Ni denotes the number of nodes with degree d\ and N 2 denote the 
number of nodes with degree d 2 , that is N\ + N 2 = N. In order to investigate the effect 
of graph structure the ratio of low and high degree nodes, i.e., N\ and N 2 is varied. 
The degrees are fixed at k\ = 5 and k 2 = 35. In Figure [2] the curves [I\ p , [I) c and [I] s 
are shown in three cases. The average degree and the standard deviation of the degree 
distribution is shown in Table |T] One can see that the new system agrees with and is 
almost indistinguishable from the compact pairwise model, in fact for bimodal graphs [7] s 
coincides with [I] c since Eqs. (l29l) - (13Ul) provide a unique solution without involving any 
approximations. Figure [7] shows that the traditional pairwise model performs relatively 
well only in the case when the standard deviation is small, that is the graph is nearly 
homogeneous. 

Consider now the case of configuration random graphs with cutoff power law degree 
distribution. These random graphs are given by a minimal degree k m i n , a maximal degree 
k m ax and a power a. The degree distribution of the graph is p{k) = Ck~° for k = 
kmin , k m i n + 1, ..., k max with the normalisation constant C given by 

kmax 

E k ~ a - 


i 

c 


k=k n 
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In Figure 0 again the functions [I] p , [J] c (t) and [/] s (t) are shown for a sparse (lower curves) 
and a dense (upper curves) power law configuration graph with power a = 2. Table [I] 
again shows the average degree and the standard deviation of the degree distribution of 
the sparse and dense networks. The value of r in both cases is r = 3'y(k)/(k 2 ). We can 
see again that the super compact PW model gives excellent agreement with the compact 
pairwise model. 


Network 

<*> 

Vtt 2 ) ~ (■ k } 2 

Bimodal 0.1 

32 

9 

Bimodal 0.5 

20 

15 

Bimodal 0.9 

8 

9 

Power law sparse 

10.1 

5.9 

Power law dense 

28.4 

26.01 


Table 1: The average degree and the standard deviation of the degree distribution of the 
graphs for which the performance of the new closure was tested. For bimodal graphs the 
degrees are k\ = 5 and = 35, the numbers in the first coloumn indicate the proportion 
of low degree nodes, i.e., N\/N. For the sparse power law graphs the degrees vary between 
k m in = 1 and kmax = 35, for the dense one k m i n = 10 and k max = 140, the power is a = 2. 


4 Discussion 


In this short paper, we derived a super compact pairwise model consisting of only 4 
equations for SIS dynamics and for heterogenous networks constructed according to the 
configuration model. This represents an improvement of going from order K, where K 
is the number of distinct degrees in the network, to order one equations, namely 4. We 
note that the closure that made the reduction possible relies on the observation that the 
distribution of susceptible nodes of degree k, which is time dependent, can be related to 
the original degree distribution of the network via a simple linear relation. We note that 
the linear relation may not be the single or unique choice, more sophisticated functional 
forms could be used based on combinatorial arguments. Moreover, the closure will not only 
encompass the first and second moment of the degree distribution but also the third. The 
new super compact model gives excellent agreement with the previously derived compact 
pairwise model. 

The accuracy of the new closure can be estimated in a semi-analytic way. The numer¬ 
ical solution of the compact PW will allow to evaluate 


E = 


S 2 -S i 

Sr 2 


Q, 


which quantifies the performance of the newly derived closure, upon using the compact 
PW model as a benchmark. Moreover, it can be shown analytically that the difference 
| [7] s (i) — [7] c (t)| can be estimated by a constant multiple of E. Further work on this model 
will include a detailed bifurcation study of the closed super compact pairwise system and 
we will aim to determine the stability of the disease free and endemic steady states. We 
will also investigate whether an explicit formula for the endemic steady state is possible. 
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If these calculations will be tractable, the stability of the disease free steady state should 
also yield Rq or at least an i?o-like quantity for the super compact PW model. 
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Figure 1: SIS epidemic propagation on a bimodal configuration random graph: simulation 
(gray thick curve), pair-wise (black dashed), compact pair-wise (black continuous), het¬ 
erogeneous pair-wise (red continuous). The two latter curves are nearly indistinguishable. 
The parameter values are N = 1000, N\ = IV 2 = N/ 2, k\ = 5, /C 2 = 35, 7 = 1 and 
r = 3 , y(k)/(k 2 ). 
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Figure 2: The curves [I\ p (dashed), [I] c (continuous) and [I] s (circles) for a bimodal 
graph with different ratios of the number of low and high degree nodes. The upper 
curves correspond to N± = O.liV, N 2 = 0.91V, the middle ones are based on N± = 0.51V, 
N 2 = 0.5 N and the lower are for N\ = 0.9N, N 2 = 0.11V. The parameter values are 
N = 1000, k\ = 5, k 2 = 35, 7 = 1 and r = 37 (k)/(k 2 ). 



Figure 3: The curves [I] p (dashed), [7] c (continuous) and [ I] s (circles) for sparse (lower 
curves) and a dense (upper curves) power law configuration graphs. The lower curves 
belong to the sparse case with k m i n = 5 and k max = 30. The upper curves belong to the 
dense case with k m i n = 10 and kmax = 140. The power is a = 2 in both cases. The 
parameter values are N = 1000, 7=1 and r = 3 r j(k)/(k 2 ). 
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